#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stdint.h>
#include <string.h>

// 大整数类型（无符号64位，支持大范围因式分解）
typedef uint64_t bigint;

/**
 * 最大公约数（欧几里得算法）：对应Python的gcd函数
 */
bigint gcd(bigint m, bigint n) {
    if (m < n) {
        bigint temp = m;
        m = n;
        n = temp;
    }
    while (m % n != 0) {
        bigint temp = m;
        m = n;
        n = temp % n;
    }
    return n;
}

/**
 * 辅助函数f(x, c, n)：x² + c mod n，对应Python的f函数
 */
bigint f(bigint x, bigint c, bigint n) {
    // 用__uint128_t避免x*x溢出（64位乘法可能超出uint64_t范围）
    return ((__uint128_t)x * x + c) % n;
}

/**
 * 因式分解主函数：对应Python的factorization函数
 * @param n: 待分解的正整数（n>=2）
 * @param fac_fun: 分解函数指针（Floyd或Brent）
 * @param factor_count: 输出参数，因子个数
 * @return: 因子数组（升序排列，动态分配，需调用者释放），失败返回NULL
 */
bigint* factorization(bigint n, void (*fac_fun)(bigint, bigint*, bigint*), int* factor_count) {
    *factor_count = 0;
    if (n < 2) {
        fprintf(stderr, "factorization: n必须大于等于2！\n");
        return NULL;
    }

    // 栈模拟Python的列表栈（最大深度32，足够分解64位整数）
    bigint stack[32];
    int stack_top = 0;
    stack[stack_top++] = n;

    // 动态数组存储因子（初始容量8，翻倍扩容）
    int factor_cap = 8;
    bigint* factors = (bigint*)malloc(factor_cap * sizeof(bigint));
    if (!factors) {
        fprintf(stderr, "factorization: 内存分配失败！\n");
        return NULL;
    }

    // 栈循环分解
    while (stack_top > 0) {
        bigint x = stack[--stack_top];

        if (x == 2) {
            // 2是素数，插入到数组开头（保持升序）
            if (*factor_count >= factor_cap) {
                factor_cap *= 2;
                bigint* temp = (bigint*)realloc(factors, factor_cap * sizeof(bigint));
                if (!temp) {
                    fprintf(stderr, "factorization: 内存分配失败！\n");
                    free(factors);
                    *factor_count = 0;
                    return NULL;
                }
                factors = temp;
            }
            // 数组元素后移，插入2到开头
            memmove(&factors[1], factors, *factor_count * sizeof(bigint));
            factors[0] = 2;
            (*factor_count)++;
            continue;
        }

        bigint p, q;
        if (x & 1) { // x是奇数，调用指定的分解函数（Floyd/Brent）
            fac_fun(x, &p, &q);
        } else { // x是偶数，分解为2和x/2
            p = 2;
            q = x >> 1;
        }

        if (p == 1) {
            // 未找到非平凡因子，q是素数，插入开头
            if (*factor_count >= factor_cap) {
                factor_cap *= 2;
                bigint* temp = (bigint*)realloc(factors, factor_cap * sizeof(bigint));
                if (!temp) {
                    fprintf(stderr, "factorization: 内存分配失败！\n");
                    free(factors);
                    *factor_count = 0;
                    return NULL;
                }
                factors = temp;
            }
            memmove(&factors[1], factors, *factor_count * sizeof(bigint));
            factors[0] = q;
            (*factor_count)++;
        } else {
            // 分解得到两个因子，入栈（栈是LIFO，先入q再入p保证分解顺序）
            stack[stack_top++] = p;
            stack[stack_top++] = q;
        }
    }

    // 调整因子数组为实际长度
    bigint* result = (bigint*)realloc(factors, *factor_count * sizeof(bigint));
    if (!result) result = factors;

    return result;
}

/**
 * Floyd查环实现的Pollard Rho算法：对应Python的floyd函数
 * @param n: 待分解的奇数（n>2）
 * @param p: 输出参数，找到的因子
 * @param q: 输出参数，n/p
 */
void floyd(bigint n, bigint* p, bigint* q) {
    bigint c = 1;
    bigint x = 2;
    bigint slow = x;
    bigint fast = x;

    while (1) {
        // slow走1步，fast走2步
        slow = f(slow, c, n);
        fast = f(f(fast, c, n), c, n);

        bigint diff = (fast > slow) ? (fast - slow) : (slow - fast); // 绝对值
        if (diff != 0) {
            bigint g = gcd(diff, n);
            if (g > 1) {
                *p = g;
                *q = n / g;
                return;
            }
        } else {
            // 查找到环且diff=0，n是质数，返回1和n
            *p = 1;
            *q = n;
            return;
        }
    }
}

/**
 * Brent查环实现的Pollard Rho算法：对应Python的brent函数（修正逻辑错误）
 * 注：Python原代码中slow = f(slow, x, n)应为f(slow, c, n)，已修正
 * @param n: 待分解的奇数（n>2）
 * @param p: 输出参数，找到的因子
 * @param q: 输出参数，n/p
 */
void brent(bigint n, bigint* p, bigint* q) {
    bigint c = 1;
    bigint x0 = 2; // 初始x值
    bigint slow = x0;
    bigint fast = f(x0, c, n);
    bigint power = 1;
    bigint steps = 1;

    while (1) {
        bigint diff = (fast > slow) ? (fast - slow) : (slow - fast); // 绝对值
        if (diff != 0) {
            bigint g = gcd(diff, n);
            if (g > 1) {
                *p = g;
                *q = n / g;
                return;
            }
        } else {
            // 查找到环且diff=0，n是质数，返回1和n
            *p = 1;
            *q = n;
            return;
        }

        // 步长等于当前power时，更新power和fast
        if (steps == power) {
            power <<= 1; // power *= 2
            steps = 0;
            fast = slow; // 重置fast到slow的位置
        }

        // slow走1步，步长+1（修正Python原代码的参数错误：x→c）
        slow = f(slow, c, n);
        steps++;
    }
}

// 打印因子数组的辅助函数
void print_factors(bigint n, bigint* factors, int factor_count) {
    printf("factorization(%llu) = [", n);
    for (int i = 0; i < factor_count; i++) {
        printf("%llu", factors[i]);
        if (i != factor_count - 1) {
            printf(", ");
        }
    }
    printf("]\n");
}

// 测试主函数
int main() {
    bigint n = 2206637; // 测试用例：2206637 = 317 × 6961
    int factor_count = 0;

    // 1. 测试Floyd查环版本
    bigint* factors_floyd = factorization(n, floyd, &factor_count);
    if (factors_floyd) {
        printf("Floyd查环版本结果：\n");
        print_factors(n, factors_floyd, factor_count);
        free(factors_floyd);
    }

    // 2. 测试Brent查环版本
    factor_count = 0;
    bigint* factors_brent = factorization(n, brent, &factor_count);
    if (factors_brent) {
        printf("Brent查环版本结果：\n");
        print_factors(n, factors_brent, factor_count);
        free(factors_brent);
    }

    // 额外测试：分解质数（应返回自身）
    bigint prime_n = 9973; // 质数
    factor_count = 0;
    bigint* factors_prime = factorization(prime_n, floyd, &factor_count);
    if (factors_prime) {
        printf("分解质数%d的结果：\n", prime_n);
        print_factors(prime_n, factors_prime, factor_count);
        free(factors_prime);
    }

    return 0;
}